
# Set up ----------------------
library(tidyverse)
library(lmtest)
library(sandwich)
library(car)
library(lavaan)
library(ri2)
library(psych)
library(margins)

# load data 
yg <- read_rds("data/yougov.rds")
lucid <- read_rds("data/lucid.rds")

# Change DV into numeric for later use
yg <- 
  yg %>% 
  mutate(
    treat_margin_num = case_when(
      treat_margin == "w1" ~ 0L,
      treat_margin == "w3" ~ 1L,
      treat_margin == "l1" ~ 2L,
      treat_margin == "l3" ~ 3L,
      TRUE ~ NA_integer_
    ),
    treat_party_num = case_when(
      treat_party == 1 ~ 0L,
      treat_party == 2 ~ 1L,
      TRUE ~ NA_integer_
    )
  )
lucid <- 
  lucid %>% 
  mutate(
    treat_margin_num = case_when(
      treat_margin == "w1" ~ 0L,
      treat_margin == "l1" ~ 1L,
      treat_margin == "l3" ~ 2L,
      treat_margin == "l5" ~ 3L,
      TRUE ~ NA_integer_
    ),
    treat_party_num = case_when(
      treat_party == "D" ~ 0L,
      treat_party == "R" ~ 1L,
      TRUE ~ NA_integer_
    ),
    treat_reminder_num = case_when(
      treat_2016remind == "control" ~ 0L,
      treat_2016remind == "reminded" ~ 1L,
      TRUE ~ NA_integer_
    ),
  )

# Function -----------------------------------
stata_robust <- function(model) {
  coeftest(model, vcov = vcovHC(model, type="HC2"))
}


# Dependent variables ------------------
# Visualize distribution
yg_dv <- yg %>% 
  dplyr::select(Q1, Q2, Q3) %>% 
  rowid_to_column() %>%
  pivot_longer(
    -rowid,
    names_to = "DV",
    values_to = "value"
  ) %>%
  mutate(
    experiment = "YouGov",
    DV = dplyr::recode(DV, 
                       Q1 = "rightful", 
                       Q2 = "legitimate", 
                       Q3 = "fair") %>% 
      factor(levels = c("fair", "legitimate", "rightful"))
  ) %>% 
  group_by(experiment, DV, value) %>% 
  summarise(n = n()) %>% 
  filter(!is.na(value))
lucid_dv <- lucid %>% 
  dplyr::select(legitimate, rightful, fair) %>% 
  rowid_to_column() %>% 
  pivot_longer(
    -rowid,
    names_to = "DV",
    values_to = "value"
  ) %>% 
  mutate(
    experiment = "Lucid",
    DV = factor(DV, levels = c("fair", "legitimate", "rightful"))
  ) %>%
  group_by(experiment, DV, value) %>% 
  summarise(n = n()) %>% 
  filter(!is.na(value)) 

dv_distribution_graph <- 
  bind_rows(yg_dv, lucid_dv) %>% 
  ggplot(aes(x = value, y = n)) + 
  geom_bar(stat = "identity") + 
  xlab(NULL) + ylab(NULL) + 
  theme_bw() +
  facet_grid(rows = vars(experiment), cols = vars(DV))

ggsave("output/figures/dv_distribution_hist.png", dv_distribution_graph,
       width = 9, height = 6)  

# correlations across DVs
cor(yg$Q1, yg$Q2, use = "complete.obs")
cor(yg$Q2, yg$Q3, use = "complete.obs")
cor(yg$Q1, yg$Q3, use = "complete.obs")

cor(lucid$legitimate, lucid$fair, use = "complete.obs")
cor(lucid$legitimate, lucid$rightful, use = "complete.obs")
cor(lucid$rightful, lucid$fair, use = "complete.obs")

# cronbach alpha 
alpha(dplyr::select(lucid, fair, legitimate, rightful))
alpha(dplyr::select(yg, Q1, Q2, Q3))


## Covariate balance -----------------------

# T tests across conditions (sample code)
t.test(lucid[lucid$treat_margin=="w1",]$is_female, 
       lucid[lucid$treat_margin=="l1",]$is_female)
t.test(lucid[lucid$treat_party=="D",]$is_female, 
       lucid[lucid$treat_party=="R",]$is_female)
t.test(lucid[lucid$treat_2016remind=="control",]$is_female, 
       lucid[lucid$treat_2016remind=="reminded",]$is_female)

# F tests
treatment_margin_yg <- 
  lm(treat_margin_num ~ is_female + age + is_not_white + is_gop, data = yg)
treatment_party_yg <- 
  lm(treat_party_num ~ is_female + age + is_not_white + is_gop, data = yg)
treatment_margin_lucid <- 
  lm(treat_margin_num ~ is_female + age + is_not_white + is_gop, data = lucid)
treatment_party_lucid <- 
  lm(treat_party_num ~ is_female + age + is_not_white + is_gop, data = lucid)
treatment_reminder_lucid <- 
  lm(treat_reminder_num ~ is_female + age + is_not_white + is_gop, data = lucid)


# Additional models -----------------------
yg_mod <- 
  lm(legitimacy ~ treat_margin + coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = filter(yg, pid7 %in% c(1,2,3,5,6,7)))
yg_baseline_with_indpd <- 
  lm(legitimacy ~ treat_margin + coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = yg)
yg_3interaction_with_indp <- 
  lm(legitimacy ~ treat_margin*is_gop*coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = yg)
yg_3interaction <- 
  lm(legitimacy ~ treat_margin*is_gop*coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = filter(yg, pid7 %in% c(1,2,3,5,6,7)))

lucid_mod <- 
  lm(legitimacy ~ treat_margin + coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
lucid_reminder_mod <- 
  lm(legitimacy ~ treat_margin + coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin + treat_2016remind, 
     data = lucid)
lucid_3interaction <- 
  lm(legitimacy ~ treat_margin*is_gop*coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)

# Table D
stata_robust(yg_baseline_with_indpd)
stata_robust(yg_mod)
stata_robust(lucid_mod)
stata_robust(lucid_reminder_mod)

# Table G
stata_robust(yg_3interaction_with_indp)
stata_robust(yg_3interaction)
stata_robust(lucid_3interaction)

# Binary change in DV 
legitimate_change_yg <- 
  lm(is_legitimate ~ treat_margin + coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = yg)
rightful_change_yg <- 
  lm(is_rightful ~ treat_margin + coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = yg)
fair_change_yg <-
  lm(is_fair ~ treat_margin + coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = yg)

legitimate_change_lucid <- 
  lm(is_legitimate ~ treat_margin + coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
rightful_change_lucid <- 
  lm(is_rightful ~ treat_margin + coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
fair_change_lucid <-
  lm(is_fair ~ treat_margin + coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)

# Table E
stata_robust(fair_change_yg)
stata_robust(rightful_change_yg)
stata_robust(legitimate_change_yg)
stata_robust(fair_change_lucid)
stata_robust(legitimate_change_lucid)
stata_robust(rightful_change_lucid)

# Average marginal effects with different baseline
yg$treat_margin <- factor(yg$treat_margin, 
                          levels = c("l1", "w1", "w3",  "l3"))
lucid$treat_margin <- factor(lucid$treat_margin, 
                             levels = c("l1", "w1", "l3", "l5"))
yg_partisans <- filter(yg, pid7 %in% c(1,2,3,5,6,7))

yg_3interaction <- 
  lm(legitimacy ~ treat_margin*is_gop*coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = yg_partisans)
lucid_3interaction <- 
  lm(legitimacy ~ treat_margin*is_gop*coparty + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)

# Table F
yg_marg_rep <- margins(yg_3interaction, data = subset(yg_partisans, is_gop == 1))
yg_marg_dem <- margins(yg_3interaction, data = subset(yg_partisans, is_gop == 0))
lucid_marg_rep <- margins(lucid_3interaction, data = subset(lucid, is_gop == 1))
lucid_marg_dem <- margins(lucid_3interaction, data = subset(lucid, is_gop == 0))


# change baseline back 
yg$treat_margin <- factor(yg$treat_margin, 
                          levels = c("w1", "w3", "l1", "l3"))
lucid$treat_margin <- factor(lucid$treat_margin, 
                             levels = c("w1", "l1", "l3", "l5"))
# 2016 reminder 
reminder_party <- 
  lm(legitimacy ~ treat_margin + coparty + treat_2016remind*is_gop +
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
reminder_4way <- 
  lm(legitimacy ~ treat_margin*coparty*is_gop*treat_2016remind + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)

# Table K
stata_robust(reminder_party)
stata_robust(reminder_4way)

# EC support 
ec_basic <- 
  lm(ec_post ~ treat_margin + treat_2016remind + ec_pre + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
ec_change_party <- 
  lm(ec_post ~ treat_margin*is_gop + treat_2016remind*is_gop + 
       ec_pre + pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
ec_coparty <- 
  lm(ec_post ~ treat_margin*coparty + 
       ec_pre + pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
ec_3interaction <- 
  lm(ec_post ~ treat_margin*is_gop*coparty  + 
       ec_pre + pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)

# Table L
stata_robust(ec_basic)
stata_robust(ec_change_party)
stata_robust(ec_coparty)
stata_robust(ec_3interaction)

# Political knowledge interaction 
lucid <- 
  lucid %>% 
  mutate(
    pol_know_bin = case_when(
      pol_know %in% c(0,1) ~ "low",
      pol_know %in% c(2,3) ~ "med",
      pol_know %in% c(4,5) ~ "hi",
      TRUE ~ NA_character_
    ) %>% 
      factor(levels = c("low", "med", "hi"))
  )

pol_know_gop <- 
  lm(legitimacy ~ treat_margin*is_gop*pol_know_bin + coparty +
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
pol_know_mod <-
  lm(legitimacy ~ treat_margin*coparty*pol_know_bin + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)

# Table H
stata_robust(pol_know_mod)
stata_robust(pol_know_gop)
pol_know_marg <- margins(pol_know_gop, data = lucid)
pol_know_marg_r <- margins(pol_know_gop, data = subset(lucid, is_gop == 1))
pol_know_marg_d <- margins(pol_know_gop, data = subset(lucid, is_gop == 0))

# Democratic support don't load so do separately
lucid <- 
  lucid %>% 
  mutate(
    dem_live_bin = case_when(
      dem_live %in% c(1,2,3,4,5,6,7) ~ "low",
      dem_live %in% c(8,9) ~ "med",
      dem_live == 10 ~ "hi",
      TRUE ~ NA_character_
    ) %>% factor(levels = c("low", "med", "hi")),
    dem_best_bin = case_when(
      dem1_best %in% c(1,2,3) ~ "low",
      dem1_best == 4 ~ "med",
      dem1_best == 5 ~ "hi",
      TRUE ~ NA_character_
    ) %>% factor(levels = c("low", "med", "hi")),
    dem_strongman_bin = case_when(
      dem4_strongman %in% c(1,2) ~ "low",
      dem4_strongman == 3 ~ "med",
      dem4_strongman %in% c(4,5) ~ "hi",
      TRUE ~ NA_character_
    ) %>% factor(levels = c("low", "med", "hi")),
  )

dem_live_mod <-
  lm(legitimacy ~ treat_margin*coparty*dem_live_bin + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
dem_best_mod <-
  lm(legitimacy ~ treat_margin*coparty*dem_best_bin + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
dem_strongman_mod <-
  lm(legitimacy ~ treat_margin*coparty*dem_strongman_bin + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)

# Table I
stata_robust(dem_live_mod)
stata_robust(dem_best_mod)
stata_robust(dem_strongman_mod)
dem_live_marg <- margins(dem_live_mod, data = lucid)
dem_best_marg <- margins(dem_best_mod, data = lucid)
dem_strongman_marg <- margins(dem_strongman_mod, data = lucid)

# Electoral sovereignty don't load so do separately
lucid <- 
  lucid %>% 
  mutate(
    ec_pre_bin = case_when(
      ec_pre == 1 ~ "hi",
      ec_pre == 2 ~ "med",
      ec_pre %in% c(3,4,5) ~ "low",
      TRUE ~ NA_character_
    ) %>% factor(levels = c("low", "med", "hi")),
    dem2_republic_bin = case_when(
      dem2_republic %in% c(1,2) ~ "low",
      dem2_republic == 3 ~ "med",
      dem2_republic %in% c(4,5) ~ "hi",
      TRUE ~ NA_character_
    ) %>% factor(levels = c("low", "med", "hi")),
    dem3_election_bin = case_when(
      dem3_election %in% c(1,2) ~ "low",
      dem3_election == 3 ~ "med",
      dem3_election %in% c(4,5) ~ "hi",
      TRUE ~ NA_character_
    ) %>% factor(levels = c("low", "med", "hi")),
  )

dem_EC_mod <-
  lm(legitimacy ~ treat_margin*coparty*ec_pre_bin + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
dem_republic_mod <-
  lm(legitimacy ~ treat_margin*coparty*dem2_republic_bin + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)
dem_election_mod <-
  lm(legitimacy ~ treat_margin*coparty*dem3_election_bin + 
       pol_interest + is_not_white + is_college + is_female + age_bin, 
     data = lucid)

# Table J
stata_robust(dem_EC_mod)
stata_robust(dem_republic_mod)
stata_robust(dem_election_mod)
dem_EC_marg <- margins(dem_EC_mod, data = lucid)
dem_republic_marg <- margins(dem_republic_mod, data = lucid)
dem_election_marg <- margins(dem_election_mod, data = lucid)


